% 参数设置
m = 4; % α粒子的质量（原子质量单位）
Z = 92; % 重核的原子序数（例如铀）
e = 1; % 基本电荷单位（这里为了简单起见，假设e=1）
k = 1; % 库仑常数（这里为了简单起见，假设k=1）

% 初始条件
v0 = linspace(1, 5, 100); % 初始速度的不同取值
theta0 = pi/2; % 初始发射角度（垂直于径向）
r0 = 10; % 初始距离

% 循环计算不同初始速度下的轨道
for i = 1:length(v0)
    % 初始化粒子的位置和速度
    r = r0;
    dr_dt = v0(i) * cos(theta0); % 初始径向速度
    theta = theta0;
    dtheta_dt = v0(i) * sin(theta0) / r0; % 初始角速度
    
    % 存储粒子的位置和速度
    r_values = [r0];
    theta_values = [theta];
    
    % 数值求解运动方程
    dt = 0.01; % 时间步长
    t = 0;
    while r > 0.01 % 防止除以零
        % 更新位置和速度
        r_old = r;
        dr_dt_old = dr_dt;
        dtheta_dt_old = dtheta_dt;
        
        % 计算新的径向速度和角速度
        dr_dt = dr_dt_old - dt * (k * Z * e^2 / (m * r^2) - r * dtheta_dt_old^2);
        dtheta_dt = dtheta_dt_old + dt * (2 * dr_dt_old * dtheta_dt_old / r);
        
        % 更新位置
        r = r_old + dt * dr_dt;
        theta = theta_old + dt * dtheta_dt;
        
        % 存储位置和角度
        r_values = [r_values, r];
        theta_values = [theta_values, theta];
        
        % 更新时间
        t = t + dt;
    end
    
    % 绘制轨道
    plot(r_values.*cos(theta_values), r_values.*sin(theta_values), 'LineWidth', 2);
    hold on;
end

% 设置图形属性
xlabel('x');
ylabel('y');
title('Particle Orbits in Repulsive Coulomb Field');
legend('v_0 = 1', 'v_0 = 2','v_0 = 5');
grid on;
axis equal;
hold off;